#!/usr/bin/perl -w
use Getopt::Long;
use Bio::DB::Sam;
use DBI;
use Bio::SeqIO;
use constant COMMIT_INTERVAL => 100_000;
use strict;

my %bases = map { $_ => 1} qw(A C G T);
my @bases = sort keys %bases;
my $bam = '/scratch/NC0001.trim_mm2_repeat_strip_contigs_SOAP.bam';
my $fasta = '/scratch/bowtie_indexes/neurospora_crassa_OR74A_7.fa';

my $reads_db = '/scratch/200M6AAXX_4.stripped.fasta';
my @regions;
my ($debug,$force);
my $outdir = '.';
my $from = 'fasta';

GetOptions(
	   'v|verbose!'  => \$debug,
	   'b|bam:s'     => \$bam,
	   'f|fasta:s'   => \$fasta,
	   'from:s'      => \$from,
	   'r|region:s'  => \@regions,
	   'db:s'        => \$reads_db,
	   'force!'      => \$force,
	   'o|out:s'     => \$outdir,
	   );
mkdir($outdir) unless -d $outdir;
# high level API
my $sam = Bio::DB::Sam->new(-bam  => $bam,
			    -fasta=> $fasta,
			    );

my $dbfile = $reads_db . ".sqlite";
my $reads_dbh;
if( ! -f $dbfile || -z $dbfile || $force ) {
    unlink($dbfile);
    $reads_dbh = DBI->connect("dbi:SQLite:dbname=$dbfile","","",
			      { RaiseError => 1, AutoCommit => 0 });
    $reads_dbh->do(<<EOF
CREATE TABLE read_info ( 
 read_id INTEGER PRIMARY KEY,
 machid  varchar(32) NOT NULL,
 seq     char(42) NOT NULL DEFAULT ''
)
EOF
		   );
    my $sth = $reads_dbh->prepare(<<EOF
INSERT INTO read_info ( machid, seq ) VALUES (?,?)
EOF
);    
    my $in = Bio::SeqIO->new(-format => 'fasta',
			     -file   => $reads_db);
    my $i = 0;
    while( my $seq = $in->next_seq ) {
	$sth->execute($seq->id, $seq->seq);
	if( ++$i % COMMIT_INTERVAL == 0) {
	    $reads_dbh->commit;
	    warn("$i rows\n");
	}
    }
    $reads_dbh->commit;
    $sth->finish;
    $reads_dbh->do(<<EOF
CREATE UNIQUE INDEX ui_machid ON read_info (machid);
EOF
    );
    $reads_dbh->commit;
} else {
    $reads_dbh = DBI->connect("dbi:SQLite:dbname=$dbfile","","",
			      { RaiseError => 1, AutoCommit => 0 });
}
my $sth = $reads_dbh->prepare(<<EOF
SELECT seq FROM read_info WHERE machid = ?
EOF
			      );
my $i = 0;
for my $region ( @regions ) {
    next unless( $region =~ /([^:]+):(\d+)(?:\.\.|-)(\d+)/ );
    my ($chrom,$start,$end) = ( $1, $2, $3);
    print "$chrom $start -> $end\n";
#my @targets    = $sam->seq_ids;
    my %counts;
    for my $a ($sam->get_features_by_location(-seq_id => $chrom,
					      -start  => $start,
					      -end    => $end) ) {
	my $seqid  = $a->seq_id;
	my $start_aln  = $a->start;
	my $end_aln    = $a->end;
	my $target = $a->target;
	my ($tid,$tstart,$tend,$tstrand) = ($target->seq_id,
					    $target->start,$target->end,
					    $target->strand);	
	my $read_seq;
	if( $from eq 'fasta' ) {
	    $read_seq = get_read($sth,$tid);
	} elsif( $from eq 'SAM' ) {
	    $read_seq = $a->target->seq->seq;
	}
	
	my $readbase_5prime = substr($read_seq,0,1);
	unless( $bases{$readbase_5prime} ) {
	    warn("base is $readbase_5prime\n");
	} else {
#	print "$seqid $start_aln..$end_aln $tid $tstart..$tend $tstrand $readbase_5prime\n";
	    my $strand = $a->strand > 0 ? 'plus' : 'minus';
	    $counts{$strand}->{$target->length}->{$readbase_5prime}++;
	}
	last if $debug && $i++ > 10;
    }

    for my $str ( keys %counts ) {	
	open(my $fh => ">$outdir/$region.$str.table" ) || die $!;
	print $fh join("\t", qw(LENGTH), @bases), "\n";
	for my $len ( sort { $a <=> $b } keys %{$counts{$str}} ) {
	    next if $len < 18 || $len > 30;
	    print $fh join("\t", $len, 
			   map { $counts{$str}->{$len}->{$_} || 0 } @bases),
	    "\n";
	}
	close($fh);
    }
    
}
sub get_read {
    my ($sth,$readid) = @_;
    $sth->execute($readid);
    my ($seq) = $sth->fetchrow_array;
    $sth->finish;
    return $seq;
}

END {
    $sth->finish if $sth;    
    $sth = undef;
    $reads_dbh->disconnect if $reads_dbh;
    $reads_dbh = undef;
}
